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Abstract. We use transport techniques to calculate the trispectrum produced in 
multiple-field inflationary models with canonical kinetic terms. Our method allows 
the time evolution of the local trispectrum parameters, tnl and <7nl, to be tracked 
throughout the inflationary phase. We illustrate our approach using examples. We 
give a simplified method to calculate the superhorizon part of the relation between 
field fluctuations on spatially flat hypersurfaces and the curvature perturbation on 
uniform density slices, (, and obtain its third-order part for the first time. We clarify 
how the 'backwards' formalism of Yokoyama et al. relates to our analysis and other 
recent work. We supply explicit formulae which enable each inflationary observable 
to be computed in any canonical model of interest, using a suitable first-order ODE 
solver. 
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1. Introduction 

Cosmological inflation predicts the generation of a primordial perturbation, (, believed 
to have seeded the temperature anisotropy of the cosmic microwave background 
( "CMB" ) and the galaxy density field. This fluctuation is sensitive to the physics that 
created it, and therefore different models of inflation typically generate perturbations 
with distinct statistical properties. These properties can be observed by measuring 
their correlation functions. We expect this approach to provide the most important 
observational constraints on an era of early-universe inflation. 

What information is encoded in these correlation functions? The two-point function 
is nearly determined by the symmetries of the background, rather than the choice of 
microphysics, although useful information may be extracted from its scale dependence. 
The higher n-point functions are much less constrained, but only the three- and four- 
point functions (the "bispectrum" and "trispectrum" ) are likely to be measured in the 
near future. Canonical single field inflation predicts a bi- and trispectrum which will 
be undetectable by present-day or near-future experiments [IHHj. But if more than one 
field is active during inflation, or noncanonical interactions are present, the three- and 
four-point functions can be measured and their properties can discriminate between 
these possibilities. 

Because of their observational relevance and constraining power, these "nongaus- 
sian" effects have received considerable attention. During inflation, each comoving k- 
mode of a light scalar field receives a perturbation when the corresponding physical 
scale crosses outside the horizon. Once outside, causality forbids any exchange between 
neighbouring regions and therefore ( must be generated by reprocessing the local fluc- 
tuations. Where only a single degree of freedom ( g is relevant, this gives [21 CO E] 

C(X) = C(X) + ^/NL(C 9 2 (X) - (O) + ^Nl4(x) + • • • , (1) 

where all quantities are evaluated at the same time, and x labels a coarse-grained 
spatial position with sub-horizon details smoothed out. This local character gives 
each correlation function a very distinctive momentum dependence. At leading-order 
the bispectrum has only one possibility, generated by the quadratic term in ([!]). Its 
amplitude is parametrized by the number /nl 0, \1Q\ , which may depend weakly on the 
smoothing scale. But the trispectrum has two possibilities, generated respectively by the 
cubic term and the square of the quadratic term. These are conventionally parametrized 
by the numbers tnl and #nl ITH - H3] . In the single-field case, tnl does not appear 
m and can be expressed in terms of /nl; the precise relation is r NL = (6/nl/5) 2 . 
Where more than one light degree of freedom is present, they may all appear in Eq. (fi]) 
and this relation is weakened to the Suyama-Yamaguchi inequality r NL ^ (6/nl/5) 2 
[T5| ITE] . The role of such relations in diagnosing the active particle spectrum during 
inflation was recently emphasized by Assassi et al. |17j . 

Transport methods. — In this paper we explain how the non-linearity parameters tnl 
and #nl can be calculated using "transport" methods. 
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Such calculations can already be carried out within the U 5N formalism" [2], [TJ, [8] , 
which requires a Taylor expansion of the background solution in small displacements 
from a chosen initial condition. An expression for tnl was given in this formalism by 
Alabidi & Lyth [13]. A comparable result for ^nl was provided by Sasaki, Valiviita & 
Wands [XT] in the context of a curvaton model, and later generalized to an arbitrary 
number of light fields in Ref. [5]. The "5iV" Taylor expansion leads to concise and 
attractive analytic results. But it is not ideally suited to numerical implementation, 
because it relies on extracting small variations which can easily be swamped by numerical 
noise. 

In Ref. [12] it was explained that the Taylor expansion can be understood as a 
variational method to compute Jacobi fields for the flow of inflationary trajectories in 
phase space. These fields can be used to explore local properties of any flow, and were 
introduced by Jacobi in his reformulation of Hamiltonian mechanics into what is now 
Hamilton-Jacobi theory. In inflation, they represent the geometrical structure which 
underlies perturbation theory in the long wavelength limit. They recur in many areas of 
physics (see, e.g., Refs. [191 ED] ) , and have been much-studied in WKB approximations 
to the path integral [2"TH2"3] . 

The Jacobi fields are the necessary ingredient to compute tnl and g-$L, but it 
is not necessary to use variational techniques to compute them. Their evolution 
can be determined equally well using an ordinary differential equation — the 'Jacobi 
equation' [21]. The equivalence was emphasized by DeWitt-Morette [21]. The Jacobi 
equation is usually preferable for numerical implementation. It can be solved using 
conventional ODE techniques and is usually much more stable against numerical 
noise. Jacobi methods are widely used in other applications, including gravitational 
lensing [251126]. 

With this motivation, one can ask whether it is possible to replace the U 5N" 
Taylor expansion with an approach based on the Jacobi equation. To do so, one gives 
an evolution equation for each n-point function. Such equations were introduced in 
Refs. [23 12B] and were originally framed in real spacejf] Real-space methods are adequate 
if one wishes to extract only the local part of the three-point function. But if one wishes 
to include more general momentum-dependence or study n-point functions for n ^ 4, 
where it is necessary to distinguish between "squeezed" and "collapsed" configurations, 
one must revert to Fourier space. In Ref. [18] it was explained how to formulate evolution 
equations for the full k-space correlation functions, which can be integrated using an 
approach similar to the "line of sight integral" used to simplify solution of the Boltzmann 
equation. In Ref. [H] this was used to give formal but explicit expressions for the n-point 
functions in terms of the Jacobi fields and their derivatives, and hence to demonstrate 
equivalence with the variational U 5N formalism" up the three-point function. 

| A similar formalism had been introduced earlier by Yokoyama, Suyama & Tanaka [29l [30] , who gave 
evolution equations for the Taylor coefficients of the SN formalism rather than the n-point functions 
directly. It was shown in Ref. |18j that these formalisms are equivalent up to the 3-point function. In 
§||3|[5|we extend this equivalence to the 4-point function. 
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In this paper we specialize this method to the trispectrum. We write a transport 
equation for the four-point function of field fluctuations 5(p a defined on spatially flat 
slices. As in Ref. [H], this can be integrated in terms of Jacobi fields and reproduces 
the variational formulae discussed above. In a second step, we express the correlation 
functions of ( in terms of those of the 5<p a . At this point the required values of r NL 
and #nl can be extracted. However, our method is not limited to obtention of the £ 
correlation functions and can be deployed to determine the correlation functions of both 
( and any isocurvature modes. 

Outline. — In £j2] we introduce the transport framework and extend it to third order. In 
§2.2| we write down the full k-dependent equation which evolves the four-point function 
on superhorizon scales. By studying the momentum-dependence of this equation, we 



can extract (in £2.3) the coefficients of the "squeezed" and "collapsed" configurations. 
We give separate evolution equations for these. 

In §[3] we demonstrate that the transport (Jacobi) method is equivalent to the 
familiar Taylor expansion of the separate universe formalism. We use our evolution 
equation to derive ordinary differential equations which evolve the separate-universe 
Taylor coefficients forward in time, and which supply the basis of an efficient numerical 
implementation. In §|4] we finish the task of extracting tnl and #nl by computing the 
relationship between ( and the field fluctuations 5ip a . Our final expressions are given 
in §4.2| We supply explicit expressions which enable each inflationary observable to be 
computed in any canonical model of interest, using a suitable first-order ODE solver. 

In §[5] we describe the alternative backwards transport method introduced by 
Yokoyama et al., and extend it to accommodate the trispectrum parameters. We 
briefly comment on the relative advantages of each formulation. In §[6] we discuss some 
representative numerical results. Finally, we conclude with a short discussion in £j7| 

Notation and conventions. — We set c = h — 1 and work in terms of the reduced 
Planck mass, Mp = (87rG) _1 where G is Newton's constant. The species of light scalar 
fields are indexed by Greek labels a, (3, . . . , . 

2. Transport Equations 

After smoothing on a length scale k^ 1 ^> (aH)^ 1 , the field value in each smoothed region 
of the universe ( "patch" ) will evolve independently, as though it were in a homogeneous 
and isotropic separate universe. Making use of the slow-roll approximation, and 
assuming that all fields are canonically normalized and minimally coupled to Einstein 
gravity, each smoothed field ip a evolves according to [TH [27J |28] 

up to gradient-suppressed corrections. In writing §2§ we have used the e-folding number 
dN = H dt as a time variable, and t is cosmic time. The index a labels the species of 
light scalar fields and u a can be interpreted as a flow vector describing the trajectory 
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of the smoothed field in phase space. In this paper we will take these indices to be 
contracted using the flat metric 5 a p, so that index placement is immaterial. 

If desired the slow-roll approximation could be abandoned by passing to a 
Hamiltonian formulation. The resulting transport equations are structurally identical, 
requiring only specification of suitable initial conditions. This method was described 
in Refs. [T8l [28] and later implemented by Dias, Frazer & Liddle [31] for the purpose 
of studying D-brane models of inflation. In this paper we will restrict ourselves to the 
slow-roll approximation, but our evolution equations are unchanged by this choice and 
can be extended immediately to the full phase space. 



2.1. Jacobi equation 

The field value varies between coarse-grained patches. Picking a fiducial patch labelled 
by the spatial position x, the field in a neighbouring patch at relative position r will be 
displaced by a small amount 5(p a , 

V9 Q (x + r) ps ip a (x.) + 5ip a (r). (3) 

At a generic position, and provided the region under consideration is not too large, 
we can expect \5(p Q \ to be small in comparison with \(p a \. With these assumptions the 
evolution of 5(p a can be obtained by making a Taylor expansion of the velocity u a in 
the neighbourhood of the fiducial trajectory. Hence, 



&8<p a (r) 1 



+ yU a p 1 slv(x)]o~¥(3(r)5(p~ / (r)5(ps(r) H . 

We now exchange r for a Fourier space description. To keep the resulting equations 
compact we employ the 'primed' DeWitt index convention introduced in Ref. [IS] . In 
this notation, a compound index such as a' includes a field label a and a momentum label 
k a , and also indicates evaluation at some common time of interest N. The summation 
convention applied to a' implies integration over momentum with measure d 3 k a /(27r) 3 , 
and summation over the species a. In this notation we find 



Ua>p>6(ppi + ^-U a 'p>y ( S(fp>5(p 7 > - (8(fp>5(py) 



dN -""-^ 2! 

+ yU a >p>y5>[8<pp>5ipy5(ps> - (5ipp'5ipy5ips') 



(5) 



Eq. ([5]) is the nonlinear Jacobi equation. We have subtracted a zero-mode, which 
amounts to discarding disconnected terms in the correlation functions. The w-matrices 
contained in ^ inherit a dependence on the fiducial region x through their dependence 
on the background fields, but the resulting connected correlation functions depend only 
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on statistical properties of the ensemble of smoothed fields. Explicitly, we find 

Ua>/3' = (27r) 3 5(k a - k j g)u aj g[9?(x)] (6) 

V/3' 7 ' = (27r) 3 5(k a - kp - k 7 )w a/ g 7 [y9(x)] (7) 
Ua'pys* = (2vr) 3 5(k a - k/3 - k 7 - k 5 )« a ^[^(x)]. (8) 

Evolution of correlation functions 

The Jacobi equation ^ summarizes evolution in the ensemble of smoothed patches. The 
w-matrices can be calculated using any suitable method, such as the long- wavelength 
limit of cosmological perturbation theory or the separate-universe approximation. 
However they are obtained, they control not only the evolution of physical field 
fluctuations but also their correlation functions. 

To show this we note that for any classical observable O not explicitly depending 
on time, the time derivative of its expectation value satisfies d(0)/diV = (dO/dN), 
provided probability is conserved|§] It also applies quantum-mechanically if O is a 
Heisenberg picture field. Transport equations for the quantum case, similar to those 
we will develop here, were given by Andrews & Hall [32] and developed by Ballentine 
& McRae [33]. The classical limit was studied by Hepp [33]. 

We define the two-point function H a >p> to satisfy 

= (<W<fy>/3')- ( 9 ) 

Recall that our index convention implies that each quantity on the right-hand side 
is evaluated at the common time of interest, N. Differentiating this expression, and 
moving the time derivative inside the expectation value as discussed above, we obtain 
an evolution equation for E^/, 

^aT = \^aT^' + ^'W/- (10) 

Use of Eq. ([5]) allows the right-hand side to be rewritten in terms of w-matrices and 
correlation functions. Working to the lowest relevant order jjj] we conclude 

dY, a ipi 
dN 

where "[^ 3p.f. ]" denotes terms containing higher-order correlation functions which 



U a i~f'Yuj'0' + U^'ySy Q ' + 3p.f. ], (11/ 



have been omitted, beginning with the three-point function. Eq. (11) will be a good 



approximation whenever these higher-order correlation functions are negligible, which 

§ Technically, the probability distribution P must vanish sufficiently rapidly on the boundary of phase 
space that u a P — > there, and therefore integration by parts inside the expectation value does not 
generate any boundary terms. 

|| Retaining higher-order contributions would reproduce the 'loop corrections' of the SN formalism; see 
Refs. [HI3S1I3B]. 
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will usually be satisfied during an epoch of quasi-exponential inflation. In that case, the 
correlation functions typically order their amplitudes in powers of if 2 [37] making the 
relative error after translation to ( of order (H/Mp) 2 C 1. A similar procedure gives 
the evolution of the three-point function. We define 

= {5ip a >5ipi3>5ipy), (12) 

and the corresponding transport equation is 

u a >\>a\'i3>y + « Q 'AV s A'/3'£y 7 ' + cyclic + [ ^ 4p.f. ], (13) 



dav/3V 



dN 

where "cyclic" denotes the two cyclic permutations of each term, and "[^ 4p.f. ]" again 
denotes terms involving higher-order correlation functions which have been discarded, 



beginning with the four-point function. As for the two-point function, Eq. (13) will be a 
good approximation whenever these are negligible in comparison with the terms which 
have been retained. 



Four-point function. — Eqs. (11) and (13) were given in Ref. [18]. In this section, for 



the first time, we give the corresponding transport equation for the four-point function. 
To do so, we must distinguish carefully between the connected and disconnected 
contributions. The disconnected contributions are always present, even in the case 
of purely Gaussian statistics, and therefore provide no new information. But if the 
perturbations develop some intrinsic nongaussianity during their evolution, this is 
encoded in the connected part of the four-point function. To obtain it we subtract 
the disconnected terms from the full four-point function, and define 

fia'P'^'S 1 = {S(p a 'S(fi/3'S(pYS(ps') ~ E^'Ey^/ — E^yE^' — E^/E^/y . (14) 

In statistical language, the four-point function {Sipa'Sipp'SipySifs') is the moment, and 
the connected part (3 a 'i3>ys> is the cumulant. 
The transport equation for fla'p^'S' is 



AN 



(ua'xPxpiys' + 3 cyclic) + UWav^AW^V*' + 11 cyclic 
+ (wa'AVV^A'yEyyE^/ + 3 cyclic) + 5p.f. ]. 



(15) 



It can be obtained by various methods, including the Gauss-Hermite cumulant 
expansion used in Ref. [27] . the method of generating functions used in Ref. [28J and 
the approach described above. 

2.3. Separation of local shapes 



The transport equations (11), (13) and (15) evolve each correlation function in its 



entirety. Although they are first order ordinary differential equations, they are not 
trivial to solve because they couple the correlation functions associated with different 
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k- and species labels^]] Indeed, the coupled system can be regarded as simply a form of 
Boltzmann hierarchy. Like the hierarchy used to compute CMB anisotropies it must be 
truncated — by discarding higher-order correlation functions — if it is to be turned into a 
practical computational tool. We will see in £j3] that it admits a similar kind of formal 
solution. But if we wish only to track the evolution of the local momentum shapes, 
then we can extract simpler "flavour" equations which do not involve the continuum of 
k-modes. These are ordinary differential equations for a finite number of variables and 
their numerical solution is straightforward. 



Eqs. (11), (13) and (15) show that (at least to this order), each correlation function 
is sourced by the correlation functions of lower order. Hence, we proceed inductively: if 
the k-dependence of the two-point function is known, then it can be used to determine 
the local k-dependence inherited by the three-point function and subsequently the four- 
point function. 

Two-point function. — Since we anticipate approximate scale-invariance, we write the 
two-point function as 

S^ = (27r) 3 5(k a + k /3 )^f, (16) 

where S a/ g has dimension of [mass] 2 but is nearly independent of k a = kp. It is this k~ 3 
dependence which will be inherited by all higher n-point functions. The possible ways 
in which this inheritance can happen correspond to the possible local ("squeezed" and 
"collapsed") momentum shapes. 

We first require a transport equation for T, a p. As described above, this is a flavour- 
only matrix, carrying indices for the species of scalar fields but not momentum labels. 



Substituting ( JX6[ ) into ( JTlj ), we conclude 

U a \£\p + UpxT^Xa- (17) 



dTj a p 
dN 



This is symbolically the same equation as the full k-space transport equation, Eq. (11), 
with primed indices exchanged for unprimed ones. 

In practice, S Q( g carries a small dependence on the fc-scale at which it is evaluated. 
This fc-dependence, typically characterized by a spectral index, can also be calculated 
by transport methods; see Dias et al. [32] ■ Recently Dias, Frazer & Liddle extended this 
method to obtain the scale-dependence of the spectral index, or "running" [31 J . 

Three-point function. — Examination of the transport equation for the three-point 
function, Eq. (13), shows that in a small time interval 5N, the change to ov^v is of the 



schematic form 5a ~ (u'a + u"££)5iV, where a prime ' applied to u indicates one of the 
field-space derivatives which generate the index structure for the w-matrices. The u'a 

% For example, the four-point function with momentum labels ki, k 2 , k 3 and k 4 couples to other 
correlation functions with momenta k! +k 2 , and so on. Had we retained loop corrections, these would 
make the hierarchy considerably more complex because each correlation function no longer couples only 
to a few other isolated k-modes, but to the whole phase space of soft superhorizon modes. Handling 
this is a computational challenge. For one approach see, eg., Ref. |38j . 
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terms generate a change 5a which is proportional to the momentum- dependence already 
carried by a. Therefore this term can reorganize the amplitudes of these shapes, but 
introduces no new types of momentum dependence. New shapes are sourced only by 
the EE terms. 



Eq. (16) shows that the product EE must generate a shape of the form k~ 3 ka 3 , 
and therefore the most general structure which can be sourced during the evolution has 
the form 

a^, y D (2*mK + k, + k,) + |^ + , (18) 

where we use the notation "D" to indicate that the three-point function contains 
this term together with others which have not been written. The matrices 

0- j a\(3'y are 

symmetric under exchange of -H- 7, but need not possess further symmetries. The 



full three-point function corresponds to the sourced contribution (18) plus an unsourced 
term appearing as its initial condition. The unsourced piece is generated by quantum 
interference effects operating around the epoch of horizon exit, and typically has a very 
complicated momentum dependence [3] . However, its amplitude is small in the canonical 
models to which we restrict attention in this paper 



After substitution of (18) into the transport equation (13), we obtain an evolution 
equation for o a ^ 7 , 

da Q |/3 7 



diV 



(19) 



Eq. (19) strictly applies only when the momenta entering the correlation function are 
not too dissimilar in magnitude. This is usually an acceptable approximation for CMB 
experiments, but a more refined analysis might be required where larger hierarchies of 
scale exist. This issue is not confined to the transport framework; it applies to results 
obtained using any method, including the familiar SN Taylor expansion. 



Four-point function. — Eqs. (17) and (19) were given in Ref. [IB]. The same analysis 



applied to the four-point function shows that, in a small time interval 8N, the change 
in the connected part of the correlation function has the schematic form 

5(5 ~ {u'(5 + «"aS + u'"'EEE)5N. (20) 

As for the three-point function, the term u ; (3 is simply a shift in the amplitude of shapes 
already present in 0. The sourced contributions are now w"aE and w^EEE. Of these, 
the EES term must generate a shape of the form k~ 3 k^ 3 k~ 3 , which can be recognized 
as a ^NL-type contribution [T3] . 

The aE term is more complex, because the momentum 5-function in u" [see Eq. (Pfl)] 



reorganizes the momenta appearing in the denominators of the three-point function ( 18 ). 
Written out explicitly, this term is 

«a'AV a A'/3' 7 ' S M'5' = ( 2lT ) 3 / d3fc A d% 5{k a - k A - k M )5(k A + kp + k 7 )5(k M + k 5 ) 

v / \ (21) 

/ Q A|^ 7 Q/3|A 7 Q 7 |/3A \ 



X M «V 7.3 I pU + 7.3 1,3 + p p I 



Transport equations for the inflationary trispectrum 



10 



plus the nontrivial permutations of a', (3', 7' and 5'. The first term in round brackets, 
~ kg 3 k~ 3 , has the form of a (?NL-type contribution. But the remaining terms involve 



A;^ , and the 8-f unctions in (21) show that kA = k a + k,5. Therefore this term generates 
a different momentum shape; it is the "collapsed" configuration, which corresponds to 
a TNL-type contribution [TJ]. It follows that the most general structure sourced by time 
evolution can be written 



ft^y, = (2nf5(k a + k, + k 7 + k s ) (gj* + 3 cyclic + ^ + 11 cyclic) , 

(22) 

where the "cyclic" pieces refer to the cyclic permutations of the preceding terms. The 
matrix g a \pys is symmetric under any exchange of 0,j,5, but has no symmetries under 
permutations involving a. The matrix T a/3 ^s is symmetric under the simultaneous 
exchanges a •f-)- /3 and 7 <H- 5, giving 12 independent elements. 



Substitution of ( 22 ) into the transport equation ( 15 ) enables us to extract individual 



evolution equations for g a \pjs and T a f3\~/6- They are 

^ = Uaxgx^s + u/3\g a \xj5 + u^xg^pxs + usxga^x 



+ UaXfi^Xl^^iiS + UaX^ax^S^n-j + UaX^Q-x^S^ /j,/3 + UaXfiu^X^ny^vS (23) 
^ — UaxTx^S + UpxTaX^S + U^a^XS + UsxT a p\~,X 

+ U^x^fiaOj^xp + Usx^npO"y\Xa- (24) 



dr a/ g| 7 5 



Note that the a-dependent source terms in the second line of ( 24 ) preserve the symmetry 
under simultaneous exchange of the index pairs (a, (3) and (7,5). We have dropped the 
initial value of cta'p'-y', even though it appears in (15) as a source term and, as a matter 
of principle, could appear in /3 a '/3'y<5' with a non-negligible coefficient. In Appendix A 



we 



show that this will usually be an acceptable approximation in models with canonically 
normalized scalar fields; the initial value of a^y remains negligible provided |r/ NL | < 1 
throughout the evolution where r < 1 is the tensor-to-scalar ratio. On the other hand, 
in non-canonical models where the initial value need not be negligible it is important 
to retain this term 



3. Equivalence to Taylor expansion method 



Eqs. (23)-(24) enable us to follow the evolution of the sourced, local-mode contributions 
to the trispectrum. As we will explain in £|4j after changing variable to ( they allow us 
to calculate the observable quantities tnl and #nl- However, they are quite different in 
appearance to the familiar expressions of the U 5N formalism" £J which take the form of 
a Taylor expansion in the initial conditions [2]. 

+ Here and below, we use the term "SN formalism" to mean a Taylor expansion in the initial conditions, 
even if the quantity being expanded is not N. 
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The connexion between these methods was explored in Ref. [18]. By formally 
integrating the transport equations, in a similar way to the "line of sight" integral used 
when solving the Boltzmann equation, it is possible to demonstrate equality with the 
U SN" expressions. In Ref. [18] this analysis was given for the two- and three-point 
functions. Here we extend it to include the four-point function. 

Integrating factor. — The "line of sight integral" naturally expresses each correlation 
function in terms of the underlying Jacobi fields. We briefly recapitulate the argument 
of Ref. |18j . Without loss of generality, we write the two-point function in the form 

^a'P' — (25) 



A suitable choice for IV j/ means it will function as an integrating factor. In writing (25) 
we have introduced a new type of primed Latin index (i', f, . . .). This has the same 
interpretation as the primed Greek indices: i' carries a flavour index i and a momentum 
label kj, which range over the same values as a and k Q . However, it indicates evaluation 



at a different time N , as follows. Substitution of Eq. (25) in (11) shows that the terms 
involving iV/3' can be removed if T is chosen to satisfy 

— U a ipil [ZO) 

Comparison with Eqs shows that T a 'i> has an interpretation as a differential 

coefficient, 

d(p a , dip a (k a ,N) ^ dcp a (N) 
Fa'i' = -x — = ^— 7] — z- x =d{k a -k i ) (27) 
dip^ d(pi(ki, N ) oipi{N ) 



Eq. (27) is sometimes described as the "Jacobi map". It has a formal solution in 



terms of a path-ordered exponential 



r aV = (27r) 3 5(k Q -k i )Pexp ( / 



N 

dN' u ai (N') ) . (28) 

N 



In this expression, P denotes the path-ordering operator which rewrites its argument 
in order of position on the trajectory: objects evaluated early on the trajectory appear 
to the right of objects evaluated later. This path-ordered exponential is related to 
the inverse of the van Vleck matrix, which is equivalent to the matrix of Jacobi fields. 



Reference to Eqs. (13) and (15) shows that, in each transport equation, this choice for T 
will absorb the terms proportional to the n-point function itself. Returning to the two- 
point function and discarding higher-order contributions, it follows that the "kernel" 
Ejy can be obtained as an integral over the source terms. It is this integral over sources 
which can be compared to the "line of sight" integral for the Boltzmann equation. 

With these choices, and working to leading order, there are no sources for the kernel 
Ej/j/. Therefore it is constant, and equal to its initial condition set at horizon crossing. 
We write this constant value iSyj/. 
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Three-point function. — When this method is applied to the three-point function, it 
transpires that the kernel is sourced. Again without loss of generality, we write 



ata'P'y' — ^a'i'^ P'j'^^'k' Ai'j'k'- 

We define uvjiy = r^ a ,u a >0/yl 'p'j'T^y and obtain 



Ai'j'k' — Ai'i'k' + 



i'j'k' 



- pN 

/ dN' Ui> m ' n >(N')S m/j/ S n ' k > + 2 cyclic 

J N 



(29) 



+ 0(# 6 ). (30) 



The integration constant Ai'j'k' is the unsourced initial condition which was neglected 
above, and the estimate 0(H 6 ) for the terms we have omitted assumes that the 
correlation functions order themselves in increasing powers of H 2 as described by 
Jarnhus & Sloth [37]. Defining 



pN pN 

Fan'? = Fa'm' / u mHlf {N') dN' = (27r) 3 5(k Q - k, - k,-)^ / u mij {N') dN', 
Jn Jn 



(31) 



(2vr) 3 <5(k Q - k, - kj)T 



we conclude 

(y a '/3'y = Fa'i'T p'jir^ki Ai'j'k' + \ F a 'i'jirp'k>Ty'i'Si'k'Sj'£/ + 2 cyclic j + loops. (32) 



It can be shown that the quantity T a ij appearing on the right-hand side of (31) is equal 
tcQ 



Oil] 



dipj(N ) dip^dif^No) 



from which it follows that Eq. (32) is equivalent to the Lyth-Rodriguez Taylor expansion 
formula for the three-point function [2]. Moreover, differentiation of (31) shows that 
T ai j satisfies the evolution equation 

drv 



dN 



fijL fill. jj. 



(34) 



Four-point function. — The analysis for the four-point function is similar. We 
introduce the integrating factor T a >i', 



fS'j'5' — Fa'i'Tpij'ry'k'Ts'e'Bi'jik'e'- 



(35) 



The kernel Bi^iyi' is given by an integral over sources, as before, which are drawn 
from the lower-order n-point functions. In this case they are the two- and three-point 
functions. Keeping only leading-order terms, we find 

Bi'j'k'i' = Bi'j'k'i' + 



dN' u i , prg ,(N , )A p .j' k '(N')S g . i/ + 11 cyclic' 

N > 

+ ( / diV' Ui' q ' r ' s i(N')S q 'j'S r 'k'S s i£i + 3 cyclic j + 0(H 8 

^ J No ' ' 



(36) 



* Direct differentiation of Eq. (28) is subtle, because of the path-ordered exponential. It is simpler to 
differentiate the Jacobi equation, Eq. (26), and then solve it using T ai as an integrating factor. 
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where we have defined U{iji k 't' = ^ii a 'U a 'p'ys , ^p'j'^'y'k'^S'e'- The integration constant 
Bi/j'k'i' is the initial value of the four-point function at time N — Nq, as for the three- 
point function. Taking the four k-modes entering the four-point function to have a 
similar time of horizon exit, and the initial time No to be around this epoch, the initial 
condition was shown in Refs. jH E] to be dominated by the correlations induced by 
decay of gravitational waves into scalar quanta. It is negligible when the amplitude of 
the four-point function is sufficiently large to be observable. On the other hand, the 
initial value of the three-point function appears in the kernel A v iyw which forms part 
of the source integral (36), and need not be entirely negligible. However, as discussed 
below Eqs. (23)-(24), and in more detail in Appendix A, its contribution to r NL or g NL 
is likely no more than 0(1) for models with acceptable |/nl|- 

To relate (36) to the expressions produced by the Taylor expansion algorithm we 
must express A v iy k i purely in terms of correlations at the initial time Ao- Combining (30 ) 
and (36) we find 



dN' Ui> p > q >(N')A p/fk ,(N') 



N 



N 



No 



dN' u i/p/q/ (N')\A p/j > k > + 



(37) 



pN' 

I dN" Upi r i s i( y N"^S r ijiS s iy 

JNq 

pN' pN' 

+ / dN" Uj' r ' s r(N")S r ip>S s 'k' + / dN" Uk'r>s'{N")S r i p iS s 'ji > 

J N J No ' 

The term involving Ap'j'k' presents no difficulties. It makes a contribution to Pa'p'^'S' of 
the form 

Pa'p^'S' 3 ^a'p'q'^p'j'^yk'^s'e'-^p'j'k'^q'e + H cyclic, (38) 

where, as above, the symbol "D" indicates that the four-point function contains this 
contribution among others. The other terms in (37) are nested integrals, and divide 
into two groups. One involves a contraction between the two w-matrices, of the form 
Ui'p'q'Up'r's'- We first focus on the other two, which involve no contraction. After 
summing over perturbations there are twenty-four such terms. Consider the specific 
choice Ui'p'q'U j' r ' s ' which appears in (37). In combination with one of the terms generated 
by simultaneously exchanging i' -h- f and k' -h- £' this generates 



pN pN 

/ dN' dN" u iW {N')u fr ' sl {N")S r , p >S s , k 'S q 

JN JNq 



(39) 



in which the integrals are no longer nested. Pairing all such terms in this way generates 
the 12 cyclic permutations of indices in (39). The corresponding contribution to the 
four-point function is 



Pa'P'-y'S' 3 ^a'p'q'^ P'r's'^-f'k'^S't'Sp'r'Ss'k'Sq'l' + H Cyclic. 



(40) 



Now focus on the contracted terms Ui> p ' q 'U p ' r ' s '. Summing over the permutations 
i' — > {j', k'} is equivalent to symmetrization over {q', r', s'}. Therefore this term can be 
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combined with the Ui> q > r i s r source in (36), giving a total contribution to the four-point 
function of the form 



Pa'fi'-y'S' 3 ^a'q'r's'^l3 / j'^Yk'^S'e / <Sq'j'Sr / k'<Ss'e' + 3 Cyclic, 

where we have defined r a /q/ r / s / to satisfy 



(41) 



N 



dN' Ui'n'r's'iN') 



No 



N 



N 



(42) 



+ (r aV / dN' u i>plql (N') / dN" u p , r , s ,(N") + [q' {r\ s'}]) . 
v Jn Jn j 

As with the previous examples of T-matrices, the momentum dependence of r a ' 5 / r / s / is 
a pure 5-function. It can be converted to a pure flavour matrix by the rule 



Fa'q'r's' = (2n) 3 S(k a - k q - k r - k s )T c 



(43) 



By explicit differentiation and back-substitution, it can be shown that this flavour matrix 
satisfies the ordinary differential equation 

dr r 



aqrs 



dN 



U a /3^/3qrs + (u a ^T /3 qr T^ s + 2 Cyclic j + Uafi^T p q T ' 7r r '$ 



(44) 



We have already seen that the lower-order Taylor coefficients T a i and T a ij are determined 



by the evolution equations (26) (with primed indices exchanged for unprimed ones) 



and (34); for an extended discussion, see Ref. [18]. These equations provide an efficient 



means to compute the U 5N coefficients" numerically. 

Returning to the four-point function, we must also include the initial condition 



Repeating the steps described above, it can be shown that 



dT 



aqr 



(45) 



(46) 



aqrs dip s (N ) d<p q (N )d<p r {N )d<p a (N o y 
Therefore we have reproduced the usual Taylor expansion formulae for the trispectrum. 



Specifically, Eq. g5J) matches (8) of Ref. 0], and Eqs. g, gOJ) and gl} match (73), 
(74) and (75) of the same reference. These expressions were later given in slightly more 
generality by Byrnes, Sasaki & Wands |14| . In the formulation given by these authors, 



Eqs. (38), (40), (41) and (45) of this paper match (36) of Ref. [H]. 



4. Transformation to the curvature perturbation 

We now have the transport equations which evolve the n-point functions of the scalar 
field perturbations during inflation, up to and including n = 4. These can be obtained 



either by solving the "shape equations", Eqs. (23)-(24), or using Eq. (44) to evolve the 



T-matrices. For the latter case, the initial conditions are T ai = 5 a i at iV = A^ , with all 
other T-matrices zero there. 
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4-1. Curvature perturbation at third order 

The scalar field fluctuations are not observable by themselves. At present we have 
observational evidence only for a single primordial fluctuation — the density fluctuation, 
which is a nonlinear and model-dependent combination of the field fluctuations. The 
appropriate combination can be deduced from the displacement 5N (measured in e-folds) 
between a fixed spatially-flat hypersurface and an adjacent uniform-density hypersurface 
with which it coincides on average. This displacement is determined by the field 
configuration on the spatially flat hypersurface. Therefore ( = 5N = 8[N(ip a )], yielding 

C = N a 5cp a + yN a p(5(p a 5(pi3 - (5<p a 5<pp)) + ^Na^Scpadcp^Sip^ - (SiPaSippSipy)) H , 

(47) 

where N a = dN/d(p a and similarly for the higher derivatives. Note that these are 
ordinary partial derivatives, with all quantities evaluated at the same time: they are 
not the nonlocal variational derivatives which appear in the Lyth-Rodriguez Taylor 



expansion. In particular, we are not using the 8N formula (47) to account for 



any time dependence of the correlation functions; this is handled by the transport 



equations. Eq. (47) is used solely to obtain the relationship between the <5<p Q and 
(. There are various other ways in which this could be obtained. Malik & Wands 
gave a comprehensive discussion [43j from the viewpoint of traditional cosmological 



perturbation theory. Another approach was used by Maldacena pQ. Eq. (47) has the 
advantage that it computes the transformation only in the superhorizon limit k/aH — > 0, 
which is all we require. 

Calculation of the derivatives N a , N a p and N a p y is tedious, although 
straightforward in principle. Ref. [18] used a raytracing method which gave the relation 
a geometrical meaning. It would be interesting to apply this technique at third order, 
but it is helpful primarily for analytic and geometric intuition rather than numerical 
optimization. Ref. [27J exploited the fact that any potential is separable for first order 
displacements to set up constants of the motion, as originally done by Garcia-Bellido 
& Wands |4TJ H3] . However, this method is relatively lengthy even for the second-order 
coefficient N a p. Here we employ a simpler alternative. 

We first focus on a single trajectory and measure the number of e-folds A" 
accumulated along it. During any period where the density decreases monotonically 
we may measure A^ as a function of p. Consider the number of e-folds AA" which elapse 
between some arbitrary point on the trajectory (the "starting point") and a nearby 
hypersurface of fixed density p c . Under the slow-roll approximation, the density at the 
starting point is simply the potential energy evaluated there. Therefore we may express 
A A" as a Taylor expansion in the difference Ap = p c — V , 

AN = N(V + Ap) - N(V) = ^Ap + i^Ap 2 + i|^Ap 3 + ■ ■ ■ . (48) 

Note that the differential coefficients are ordinary derivatives taken along the trajectory. 



In Eq. (48) they are evaluated at the starting point. 
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We now perturb the starting point by an amount 5<p a while keeping the final 
hypersurface fixed. In general 5(p a will not be aligned with the inflationary trajectory 
used to construct the p-derivatives in Eq. (48 ), which therefore vary. The same is true for 
the displacement Ap. Accounting for both these effects changes the total elapsed e-folds 
by an amount 8 (AN). Finally, to study fluctuations around the hypersurface p = p c we 
take the limit Ap — > 0, after which 5 (AN) — > £. The advantage of this method is that 
it uses the handful of low-order derivatives appearing in Eq. (48) to isolate the limited 
information we require regarding local properties of the transformation: higher-order 
information is discarded at the outset. This contrasts with the constants-of- motion 
approach used in Ref. [27], where high-order information is implicitly kept through the 
majority of the computation, although it is never used. 

Under a shift of the starting point we conclude 



S(Ap) = -V a 8<p c 



2! 



(49) 



By retaining contributions to p from the kinetic energy, and evaluating the differential 
coefficients in (48) without use of the slow-roll approximation, this approach could be 
extended to provide the transformation from the full phase space variables 5<p a , 5<p a to 



£. This was done in Ref. 

Invoking the slow- roll approximation, we may calculate the derivative dN / dp, 



dN 
dp" 



dN dt 
~dt d~V 



3H 2 



1 V 



(50) 



where dV/dt is to be computed along the trajectory 7. Higher derivatives can be 
obtained in the same way, by repeated differentiation with respect to t and use of the 
chain rule to convert these into derivatives with respect to p. We obtain 



d 2 iV 



dp 2 
d 3 7V 



Mp 



- 2- 



yv a v p v a p 



dp 3 



V X V X 
(V X V X )3 



(V X V> 



A 



(V X V X 



(VxV x 



(VxV x 



N n = — 



The first and second-order variations are 

diV 
dp" 

dN d 2 N 
N a p = - — V a 



Va 



dp 



P 



dp 2 



V a V p + (V a Ap + V p A a ) , 



(51) 

)• 

(52) 

(53) 
(54) 



which agree with existing expressions in the literature. (See below for the definition of 



A a .) At third order we find 
d 3 iV 



N 



dp 3 
1 



V n V R K 



7 



d 2 N 



dp 2 



VaVfa + cyclic - 



dA^ 
dp" 



V, 



ap'y 



+ 772 (A*Vs 7 + cyclic) + 



Mp 



^2 (B a pV 7 + cyclic) + ^ (CaVpVy + cyclic) 



(55) 
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The tensors A a , B a R and C a have been defined to satisfy 



A, 



a/3 



Vg 

V X V X 

Vol 

v x v x 

Vgfi 

v x v x 



VVkV^ 

(v x v x y 

V K V Ka Vf3 + V K V K pV a , VV K V Ka V e V ef s VV Ka V R p vv K v Kafi 

— ^ , rr; r o , , „ L , % „ — Z- 



(V X V X ) 



(v x v x r 



(v x v x ; 



+ 2 



VeMnVneVe 

(y x v x y 



- 12 



VV K V Ke V e V p V pa 



(v x v x y 



+ 4 



{V x V x f 



+ 2 



(v x v x y 
vv K v Kea v e 



(v x v x y 



(56) 
(57) 
(58) 



4-2. Inflationary observables 

Finally, we must assemble all these contributions to obtain expressions for tnl and <?nl- 
We find 



P c = N a NpE af} 



:./: 



NL 



tnl 



54 
25 



0NL 



NaNpN^Nsga^ys + SNaNpN^N^a^apE^ + N a p 7 NsN x N p 



(59) 
(60) 

(61) 

(62) 



5. Alternative approaches 

In this paper, our approach to calculating the statistics of the curvature perturbation has 



been to develop transport equations for objects such as the n-point functions [Eqs. (11 ) 



(13) and (15)], or their shape tensors [Eqs. (17), (19) and (23)-(24)]. The results of ^3 
show that this is equivalent to the Lyth-Rodriguez Taylor expansion 



6cp a = T ai SLpi + —YaijdLpidLpj + — T aijk 5cpi5(pj5(p k H 



(63) 



where we recall that objects with Greek indices are evaluated at time N, and those with 
Latin indices at some earlier time iVo, which is usually taken as the common time of 
horizon exit for the k-modes under consideration^] 

'Forward' and 'backward' methods. — To solve these equations we must supply a 
boundary condition at iV = N , but we are free to choose how this is done: we may start 
either with Nq at the initial epoch and evolve N forward to the time of interest, or fix iV 
at this time and evolve N backwards. These approaches are distinct but equally valid, 
because (at least during inflation) there is no obstacle to computing the relevant initial 
conditions at any time of our choosing. The transport equations we have described in 
this paper are of the forwards variety. 



Indeed, one can verify that inserting (63 1 into i5l and equating coefficients order-by-order reproduces 



the F-matrix evolution equations with unprimed indices [Eqs. (26), (34) and (44)] 
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Eq. (63 ) shows explicitly what must be computed in order to completely characterize 
the fluctuations 5(p a at any given order. At first order in a <i-field slow-roll model, we 
require the d 2 independent components of the Jacobi map T a i. At second order there are 
d 3 components of T ai j, reduced to d 2 (d + l)/2 after accounting for symmetries. Finally 
at third order there are d 4 components of T ai j k , which reduce to d 2 (d + l)(d + 2)/6 
independent components after symmetries. We conclude that to compute all two-point 
functions in such a model requires solution of 0(d 2 ) differential equations. Likewise, 
all three-point functions requires 0(d 3 ) equations, and all four-point functions requires 
0(<i 4 ) equations. This is to be expected, because there are 0(d m ) independent m-point 
functions. 

Autocorrelation functions of £ only. — Sometimes we do not require all correlation 
functions, but only the autocorrelation functions of £. In such cases it would be 
advantageous if an autonomous set of transport equations could be set up for the Taylor 
coefficients of ( rather than 5(f a , 

5N = Nidtpi + ^Nijdcfidtpj + ^%fc<fyi^-<fy>fc. (64) 

This would require the solution of only 0(<i m ~ 1 ) independent equations to obtain the 
m-point function of (. This saving could be helpful in models with a large number of 
fields. There is currently no forwards formulation of this type but a set of backwards 
equations were given by Yokoyama, Suyama & Tanaka [29J [30], and later extended to 
the trispectrum [4"5] . 

The Taylor coefficients for N can be expressed in terms of the T-matrices, 

Nt = N a T ai (65) 

Nij = N a Taij + NapTaiTpj (66) 

Nijk = N a T a ijk + N a /3^T + (NapFpiTajk + 2 cyclic) . (67) 

We could attempt to obtain forward transport equations by direct differentiation with 
respect to time followed by use of the T-matrix evolution equations. But this does 
not generate a closed set of autonomous equations because derivatives of the iV a ... also 
appear, which obstruct an attempt to eliminate the T-matrices in favour of their N 
counterparts. 

Instead, the backwards equations of Yokoyama et al. can be derived as follows. As 
described above, we fix N to be the late time of interest and aim to evolve No backwards. 



The backwards evolution of T ai can be obtained very simply by differentiating (63 ) while 



keeping 5tp a fixed, or alternatively by differentiating (28) with respect to A^ . Whichever 



method is chosen, we find dr^/diVo = — T a jUji. Subsequently differentiating (65) with 



respect to N and using this relation, we obtain an autonomous set of equations for Aj. 



(68) 
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This technique can be extended to higher orders, giving evolution equations for Nij and 
A^fff] We find 

= -N k u kij - N ik u kj - N jk u ki , (71) 

and 

-N £ u ti jk ~ (Nau ejk + N jke u u + 2 cyclic) . (72) 



dN ijk 



dN 

The first of these was given in Ref. [T5]. Here we have extended the method to include 
Nijk, which enables trispectrum quantities to be calculated. These equations should be 
solved with initial conditions chosen so that iVj, iVy and N,ij k equal the transformation 
matrices N a , N a p and N a ^, respectively, at N = N. 

If we require only the bispectrum of ( and are prepared to take the field fluctuations 
at time iVo to be Gaussian and uncorrelated with each other, then more is possible. 
Under these circumstances, Yokoyama et al. showed that the 0(d 2 ) equations for 
could be replaced by only 0(d) equations for an auxiliary quantity O a = T ai Ni [2"9"| IBP]. 

Constraint for first-order coefficients. — There is a further simplification which can 
be made for the Ni system. Using the flow equation d<pi = Ui dN, it follows that the 
displacement d<pi = precisely tangent to the trajectory generates a change in the 
e-foldings required to reach the final uniform density slice corresponding to 

5N = Ui Ni = -1. (73) 

This implies that one of the iVj can be determined algebraically in terms of the others, 
without solving a separate differential equation. Therefore, in a two-field model, the 



Yokoyama et al. equations (68) can be decoupled 



iN *-f n \ x ,-u H ]N,+ n ^, (74) 



where we have labelled the fields and x- A similar equation can be given for N x , but 



it is unnecessary because (73) can be used to obtain N x once is known. Although 
the possibility of decoupling these equations is interesting, it confers no particular 
advantages. 

A variation of the Yokoyama et al. formulation was recently given by Mazumdar & 
Wang [IB] in which they pointed out the possibility of this decoupling in the two-field 



case, although without making explicit use of the constraint (73). Their analysis is 



equivalent to the one presented here, and in Appendix A of Ref. (XH3 - Mazumdar & 
ff The evolution equations for T ai j and T ai j k , which are required to obtain these results, are 

" am Uj mtj aim Uj mj x amj ^mi V u;7 / 



dr Q 



dN 



^ am^mijk (^aim^rnjk H~ -T amj k^mi 2 Cyclic^ . (^0) 
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Wang ascribed the possibility of decoupling to the choice of coordinates used in their 
derivation. However, the evolution equation (68) can be derived using any convenient 
method and is independent of such choices. The argument above shows that decoupling 
is a consequence of the constraint (73), and is a special feature of the two-field system. 
In a general <i-field model, the best that can be obtained is a coupled system of d — 1 
equations. 



6. Numerical results 



We now illustrate the transport approach using a number of concrete models. For each 
model, we numerically solve Eqs. (17), (19) and (23)-(24), and use Eqs. (60)-(62) to 
determine the values of /nl, tnl and #nl from horizon crossing onwards. We label the 
number of e-folds of inflation from N = at horizon exit. 



6.1. Numerical Examples 

D-brane model. — Our first example was studied by Dias, Frazer & Liddle |31j . 
It is an approximation to inflation driven by the motion of a D-brane in a warped 
throat, allowing for angular degrees of freedom. In that study, the authors employed 
the transport approach to calculate the distribution of observable parameters over a 
large number of realizations of their model. However, they restricted attention to the 
spectrum and local-type bispectrum. Here we present the evolution of the local-type 
trispectrum parameters for one typical realization. 
The potential is given by 

V = a + «i0i + a 3 (f)l + (3(f)2, (75) 

which contains an inflexion point in the 0i direction. Inflation occurs close to this 
inflexion point. We choose a = 100M 2 M|, a x = M 2 M P , a 3 = 5M 2 /M P , f} = 5M 2 M P , 
^lexit = 0.5M P , and </>2exit = 0.5M P , where the subscript 'exit' indicates these are the 
initial values of the fields at horizon exit. M is an overall normalisation, which can 
be fixed to match the WMAP normalization of the power spectrum. These initial 
conditions have been chosen to give 60-efolds of inflation, taking inflation to end when 
e = 1. Allowing the system to evolve past this point would lead to erroneous results 
because we are employing slow-roll equations of motion. As explained in £j2]this could be 
resolved by writing transport equations in the full phase-space. However, for simplicity, 
we do not do so here. In Fig. [l]we give the evolution of tnl, ^nl, and (6/nl/5) 2 for this 
choice of parameters and initial conditions. 

For single-field models we recall that tnl = (6/nl/5) 2 , which is relaxed to an 
inequality in multiple-field models [15, 16J. The use of the relative magnitude of r NL and 
(6/nl/5) 2 as a diagnostic of the spectrum of active fields during inflation was emphasized 
by Smidt et al. [J7], who made a forecast of observational prospects. Very recently, 
Assassi et al. [T7| gave precise formulae in terms of the spectrum of single-particle 
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states. This signature of multiple active fields is clearly visible in Fig. [TJ although in 
this realization the nongaussian parameters are too small to be observable. (As a point 
of principle an inflexion point potential may give rise to a large local bispectrum [IE] 
and trispectrum |l9] , via the hilltop mechanism suggested by Kim et al. [50] . However, 
an observable signal can usually be obtained only for finely tuned initial conditions and 
parameter choices.) 
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Figure 1. Evolution of tnl, 3nl and (6/5/nl) 2 for the inflexion-point potential (75 1 
Initial conditions and parameter choices are described in the main text. 



Quadratic-exponential model. — Our second example was constructed by Byrnes 
et al. [51] as an example of a product-separable model which could give rise to a large 
/nl for finely-tuned initial conditions. It was later studied by Elliston et al. [52] and 
Huston et al. [53] . 
The potential is 

V = M^je-^l (76) 

We choose the parameter values and initial conditions A = 0.05/M|, 0i ex it = 16M P , and 
02exit = O.OOIMp, and fix M as before to match the WMAP normalization. These initial 
values also give 60-efolds of inflation. They have been chosen to select a background 
trajectory which gives rise to significant nongaussianity. In Fig. [2j we present the 
evolution of the tnl and ^nl parameters in this model for the first time. We also show 
the evolution of (6/nl/5) 2 . However, although the /nl and r NL parameters are large 
at the end of inflation, it is important to note that the fluctuations are still evolving 
at this time. Therefore the model is not predictive by itself: it must be supplemented 
by post-inflationary evolution, which tracks the fluctuations until the surface of last 
scattering, or explains how all isocurvature modes eventually decay. 
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Figure 2. Evolution of t nl , <? nl and (6/5/nl) 2 for the potential (76). Initial 
conditions and parameter choices are described in the main text. 



Non-separable hybrid model. — Finally, we present results for a hybrid-type 
potential in the large field regime studied by Mulryne, Orani & Rajantie [54J. This is an 
example of a non-separable potential. For general initial conditions, no analytic estimate 
is known for any of /nl, tnl or <7nl, even assuming slow-roll. Therefore numerical 
methods, such as our implementation of the transport equations, become essential. The 
potential contains a hilltop region, and parameter choices and initial conditions can be 
chosen so that the model is of the type discussed by Kim et al. [50J. This gives rise to 
large nongaussianity for initial conditions sufficiently close to the hilltop. 
The potential satisfies 



V = M 



(77) 



and we choose the parameter values g 2 = t> 2 /0 2 rit , m 2 = v 2 , v = 0.2M P , crit = 20M P 
and A = 5. The initial conditions are 0i ex it = 15.5Mp and 02exit = 0.0015Mp. 
As above, these initial values give 60-efolds of inflation and have been adjusted to 
produce significant nongaussianity. M is adjusted as before. In Fig. |3j we present the 
evolution of the tnl and #nl parameters in this model for the first time. In contrast 
to the previous example, the statistics here approach constant values before the end of 
inflation, reflecting the fact that isocurvature modes decay. We also give the evolution 
of (6/ NL /5) 2 . 

7. Discussion and Conclusions 



In this paper we have provided transport equations to evolve the four-point functions of 
a collection of light scalar fields during an inflationary phase. The transport system can 
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Figure 3. Evolution of t nl , <? nl and (6/5/nl) 2 for the potential (77). Initial 
conditions and parameter choices are described in the main text. 



be thought of as a form of Boltzmann hierarchy, and can be solved by similar methods. 
Since inflationary fluctuations are typically close to Gaussian, connected correlation 
functions of increasing order are typically decreasing in amplitude. Therefore only a 
few low-order functions are important in sourcing those of higher order. Truncating 
the hierarchy to include only these sources generates the local-type "squeezed" and 
"collapsed" configurations. We parametrize the amplitude of these configurations 
with "shape tensors" for which we have supplied evolution equations. Expressing the 
correlation functions of ( in terms of those of the 5<p a , it is possible to extract tnl and 
<7nl- This analysis was given in §[2] 

This method of integrating the transport hierarchy expresses the correlation 
functions in terms of the Jacobi fields generated by the underlying phase space flow, 
and their derivatives. One can regard this as a statement of the separate universe 
approximation. The "Jacobi map" relates these fields to the variation of a general 
solution of the equations of motion with respect to its constants of integration. Using 
this equivalence, we have shown that the result reproduces the familiar Taylor expansion 
used by Lyth & Rodriguez. The procedure can be viewed as an application of classical 
Hamilton-Jacobi theory. 

Our equations supply a toolkit which can be used to study the evolution of 
inflationary observables in any multi- field model of interest, provided all fields possess 
canonical kinetic terms. There are two equivalent approaches. First, one can solve 



Eqs. (17), (19), (23) and (24) for the shape tensors corresponding to the two-, three- and 



four-point functions, using suitable initial conditions. Eqs. (59)-(61) can then be used 



to extract observables. Alternatively, one can solve the evolution equations (26) (after 



exchanging primed for unprimed indices), (34) and (44) for the Taylor coefficients of the 
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U 5N formalism", applied to the field fluctuations. Once these are known, Eqs. (65)-(67) 
can be used to exchange them for the Taylor coefficients of N itself. The usual formulae 
then allow observables to be computed. If the spectral index or its running are required, 
they can be extracted using the methods described by Dias et al. [SU [39] . 

Assuming slow-roll, either method requires the solution of 0(d m ) equations to 
obtain the m-point functions of a <i-field model. Since there are 0(d m ) independent 
correlation functions it will not be possible to reduce this asymptotic complexity. But if 
only the autocorrelation functions of ( are required, then it may be advantageous to use 
the 'backwards' formalism introduced by Yokoyama, Suyama & Tanaka, in which one 
can reduce the number of equations to be solved to Old" 1-1 ) by forfeiting the possibility 
of obtaining correlation functions with insertions of isocurvature modes. [For clarity, we 
emphasize that the formalism of Yokoyama et al. correctly accounts for the influence 
of these isocurvature modes on the evolution of the £ correlation functions. But it is 
not possible to determine mixed correlation functions, such as ((s), where s is a field 
space direction orthogonal to £.] Unfortunately, it is often necessary to know something 
about such correlation functions to determine whether unquenched isocurvature modes 
remain, which could change the inflationary prediction by transferring their energy to 
the curvature fluctuation during or after reheating. (We refer to Ref. [18] for a more 
comprehensive discussion.) But in some cases this may not be a concern, and where this 
is true our extension of the formalism of Yokoyama et al. allows trispectrum parameters 
to be obtained. 
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Appendix A. Contributions to the four-point function from the initial 
condition of the three-point function 



In this appendix we verify the claim made in £2^3 that the arbitrary initial condition 
for the three-point function makes a negligible contribution to the sourced component 
of the four-point function. In the text, this was used to conclude that the initial value 



need not be retained in Eqs. (23)-(24). 



It was first proved by Lyth & Zaballa that the initial condition for the three-point 
function could be neglected in comparison with the sourced contribution whenever the 
sum of the two was large enough to be observed [ID]. Their argument was later simplified 
by Vernizzi & Wands [41] . The same result for the four-point function follows from the 
analysis of Refs. [U |6j. However, we are unaware of a similar demonstration for the 
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question addressed in this appendix — the contribution of initial value of the three-po'mt 
function to the sourced component of the four-point function. 

We work with the variational formulation of the separate universe approximation, 
as discussed by Lyth & Rodriguez [2]. We write 

C = NiStpi + -NijSpiSpj + ■■■ , (A.l) 

where the Latin indices i, j, . . . , have the same meaning as in the main text. We define 
the trispectrum to be the four-point function with its momentum-conservation 8- 
function stripped away, 

(C(k 1 )C(k 2 )C(k 3 )C(k 4 )> = (2n) 3 5(k 1 + k 2 + k 3 + k 4 )T c . (A.2) 

Using the initial value of the three-point function computed in Ref. [3j , the corresponding 
contribution to the sourced part of can be written 

tt4 ■ s U2L.2 u3 

+[k 1 ->{k 2 ,k 3 }] 

+ cyclic permutations k 4 — > {ki, k 2 , k 3 } 



(A.3) 



where "*" denotes evaluation at horizon exit, kt — ki + & 2 + fci 4 , the summation is over 
all simultaneous permutations of the index set {(3, 7, e} and the momenta {k 2 , k 3 , k 14 }, 
and we have defined kn = | ki + k 4 1 . 

This contribution can be divided into an effective g^L, an effective tnl, and an 
'equilateral-type' term which does not fit naturally into either of the local-type shapes. 
The effective #nl can be written 

(the placement of indices is immaterial in this and other expressions, since contraction 
occurs under the Kronecker-5) , and the effective r NL is 

1 NiNijNj tk . 

These expressions can be simplified. Introducing the scalar-to-tensor ratio r and the 
spectral index n s , we find 

25 

A#nl = TT^r(n s - 1 + 2e*) < 1 (A.6) 
1152 

Ar NL = -— r/NL, (A.7) 
6b 
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where /nl is the sourced local-mode contribution to the three-point function. The (?nl 
contribution is clearly negligible. The tnl contribution is negligible provided |r/ NL | < 1. 
Taking the bound on r to be roughly r < 0.1, this term can be observationally relevant 
only if |/nl| ^ 50. This is already on the verge of being ruled out by experiment, so the 
r NL contribution is likely to be no more than 0(1) in most acceptable models. It could 
perhaps be kept if very accurate estimates are required. 
Finally, the equilateral-type term is 



„2p 

"2^3 . //„ 2 , ;2^ 



i' 8 ^ 2 + ^ ~ ^ + h) ~ k2h{k2 + k 



3; 



+ [kx -> {k 2 , k 3 }] + (cyclic k 4 -> {k 1; k 2 , k 3 }) 

(A.8) 



The coefficients of these contributions are related to those of Eqs. (A. 4) and ( A.5[ ), and 
will therefore not typically be large. 
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